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£ . ABSTRACT 

d \ With the advent of high-quality surveys in cosmology the full three-point correlation function 

will be a valuable statistic for describing structure formation models. It contains information 
on cosmological parameters and detailed halo properties that cannot be extracted from the 
CNI | two-point correlation function. We use the halo clustering model to analytically calculate the 

three-point correlation function (3PCF) for general cosmological fields. We present detailed 
results for the configuration dependence of the 3-dimensional mass and galaxy distributions 
and the 2-dimensional cosmic shear field. We work in real space, where higher order correla- 
tion functions on small scales are easier to measure and interpret, but halo model calculations 
get rapidly intractable. Hence we develop techniques that allow for the accurate evaluation 
. of all the contributing terms to real space correlations. We apply them to the 3PCF to show 

■ how its configuration and scale dependence changes as one transits from the nonlinear to the 
psj | quasilinear regime, and how this depends on the relative contributions from the 1-, 2- and 
f^) . 3 -halo terms. 

The 3PCF violates the hierarchical ansatz in both its scale and configuration depen- 
dence. We study the behavior of the coefficient Q in the expansion: £( r i2j r 23, ?*3i) = 
Q[C( r u)C( r 23) + £( r i2)£(?~3i) + £(?*23)f (V31)] from large, quasilinear scales down to about 
20 kpc. We find that the nonlinear 3PCF is sensitive to the halo profile of massive halos, 
especially its inner slope. We model the distribution of galaxies in halos and show that the 
3PCF of red galaxies has a weaker configuration and scale dependence than the mass, while 
for blue galaxies it is very sensitive to the parameters of the galaxy formation model. The 
3PCF from weak lensing on the other hand shows different scalings due to projection effects 
and a sensitivity to cosmological parameters. We discuss how our results can be applied to 
rS | various analytical calculations: covariances of two-point correlation function, the pairwise 

■ peculiar velocity dispersion, higher-order shear correlations, and to extend the halo model by 

including the effects of halo triaxiality and substructure on statistical measures. 

Key words: cosmology: theory — large-scale structure of universe 
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1 INTRODUCTION 

Understanding the nature of the large-scale structure and the evolutionary history of the universe is the central aim of cosmology. Therefore, 
quantifying the observed distribution of matter rigorously is of fundamental importance, and n-point correlation functions have been the 
most widely used statistical tools for this purpose (e.g., see Peebles 1980). The inflationary scenario predicts Gaussian initial conditions. The 
statistical properties of the Gaussian field are fully characterized by the two-point correlation function (2PCF) or its Fourier transform, the 
power spectrum. Based on this idea, the two-point statistical measures for various cosmological fields have been extensively used for testing 
the paradigm as well as constraining cosmological parameters; the galaxy distribution in the redshift catalog (e.g., Totsuji & Kihara 1969; 
Hamilton & Tegmark 2002), the cosmic microwave background anisotropy (e.g., de Bernardis 2000) and the cosmic shear field (e.g., Van 
Waerbeke et al. 2001). However, even for Gaussian initial conditions, nonlinear gravitational instability induces non-Gaussian signatures in 
the mass distribution, which contain information on the nature of gravity and the dark matter. The three-point correlation function (3PCF), or 
equivalently, its Fourier-transformed counterpart, the bispectrum, is the lowest order statistical tool to probe the non-Gaussianity. The 3PCF 
can place strong constraints on models of structure formation. As an example, the measurements of the bispectrum in the galaxy distribution 
at the large scales ( ^ 10 h~ Mpc) have been used to determine the bias parameter as well as to constrain the primordial non-Gaussianity 
(e.g., Feldmann et al. 2001; Verde et al. 2001). It is also expected that the 3PCF of the cosmic shear field can be used to precisely determine 



* E-mail: mtakada@hep.upenn.edu 
f E-mail: bjain@physics.upenn.edu 



2 M. Takada & B. Jain 



the cosmological parameters complementary to the measurements of the 2PCF (e.g., Bernardeau, Mellier & Van Waerbeke 2002; Bernardeau, 
Van Waerbeke, & Mellier 2003). 

In this paper, we focus on developing a theoretical model of the 3PCF in real space for 3D and 2D cosmological fields. In practice, on 
small scales the 3PCF would be easier to measure from observational data over the bispectrum, since it does not require the Fourier transform 
for the survey data that usually have a complicated geometry of data fields. 

Theoretical models of the weakly nonlinear 3PCF have been well studied in the literature based on perturbation theory (Fry 1984b; 
Jing & Borner 1997; Gaztafiaga & Bernardeau 1998; Frieman & Gaztanaga 1999; Barriga & Gaztanaga 2002). Perturbation theory can 
describe properties of the dark matter and galaxy clustering on large scales £ 10 h _1 Mpc and predict that the 3PCF depends on the shape 
of triangle configuration and, as a result, contains information of the primordial power spectrum and the galaxy biasing. Historically, the 
pioneering measurement of the galaxy 3PCF done by Peebles & Groth (1975) (also Groth & Peebles 1977) proposed the "hierarchical form", 
C(fi2, »"23, r 3 i) = Q[£(ri2)£(r23) + C( r 23)C( r 3i) + £( r 3i)£( r i2)], with the constant Q ~ 1.3. However, subsequent work has revealed 
that the measured 3PCF does not obey the hierarchical form rigorously and the large-scale amplitudes can be explained by the perturbation 
theory results of the cold dark matter (CDM) models, if the biasing relation is correctly taken into account for the analysis (Jing & Borner 
1998; Frieman & Gaztanaga 1999). 

On the other hand, a quantitative theoretical model of the 3PCF in the strongly nonlinear regime is still lacking except for studies relying 
on TV-body simulations (Matsubara & Suto 1994; Suto & Matsubara 1994; Jing & Borner 1998) 1 . Simulations provide only limited physical 
insight into the complex non-linear phenomena involved in gravitational clustering, and are intractable for performing multiple evaluations 
in parameter space. 

Therefore, the main purpose of this paper is to develop an analytical model for predicting the 3PCF applicable to both the linear and 
nonlinear regimes. For this purpose, we need a model to correctly describe the redshift evolution and statistical properties of gravitational 
clustering up to the three-point level. We employ the so-called dark matter halo model, where gravitational clustering is described in terms 
of correlations between and within dark matter halos. Originally, this model was developed to express nonlinear clustering as the real-space 
convolution of halo density profiles (Neymann & Scott 1952; Peebles 1974; McClelland & Silk 1977; and also see Scherrer & Bertschinger 
1991; Sheth & Jain 1997; Yano & Gouda 1999; Ma & Fry 2000b,c). Most recent works have relied on the the Fourier-space formulation, 
since the forms of the power spectrum and the bispectrum become much simpler (Seljak 2000; Ma & Fry 2000c; Peacock & Smith 2000; 
Scoccimarro et al. 2001; Cooray & Hu 2001a,b; Berlind & Weinberg 2002; Hamana, Yoshida & Suto 2002; Takada & Jain 2002 hereafter 
TJ02; Scranton 2002; and also Cooray & Sheth 2002 for a recent review). The halo model appears remarkably successful in that, even though 
it relies on rather simplified assumptions, it has reproduced results from numerical simulations (Seljak 2000; Ma & Fry 2000c; Scoccimarro 
et al. 2001; TJ02) and also allowed for interpretations of observational results of galaxy clustering (Seljak 2000; Scoccimarro et al. 2001). 

We formulate the 3PCF model so that it can be applied to general 3D and 2D cosmological fields, such as the mass and galaxy distri- 
butions and the cosmic shear fields. Our method is built on the real-space formulation for the correlations of three particles in one halo. This 
is because the real-space approach enables us to compute the one-halo contribution to the 3PCF by a 4-dimensional integration, which is an 
advantage compared to the Fourier space approach. For the 2- and 3-halo terms, we rely on the Fourier-space approach and the approxima- 
tions presented in Scoccimarro et al. (2001; also see TJ02). We study the transition from the quasi-linear to nonlinear regimes and the relative 
contribution of the different terms to the 3PCF. We show the halo model predictions for the 3PCF of the mass and galaxy distribution and 
the weak lensing convergence field for the currently favored CDM model. To do this, we will focus on the dependences of the 3PCF on the 
triangle configurations as well as on the properties of the halo profile, its inner slope and concentration. 

This paper is organized as follows. In §2 we present the formalism to compute the 3PCF of the 3D density field based on the real-space 
halo model and review the Fourier-space halo approach and halo model ingredients (halo density profile, mass function and bias model) we 
will use in this paper. In §3, we present the method to compute the angular 3PCF. In §4 we show the results of the halo model predictions for 
the 3PCF of the mass and galaxy distributions and the weak lensing convergence field. Finally, §5 is devoted to a summary and discussion. 
We give some useful approximations for calculating the 3PCF in Appendix A and B. We use the currently favored CDM model (ACDM) 
with Q m o = 0.3, Q,\o = 0.7, h = 0.7 and as = 0.9. Here f2 m0 and Qxo are the present-day density parameters of matter and cosmological 
constant, h is the Hubble parameter, and as is the rms mass fluctuation in a sphere of 8h _1 Mpc radius. The choice of as for this model is 
motivated by the cluster abundance analysis (Eke, Cole & Frenk 1996). 



2 DARK MATTER HALO APPROACH 
2.1 Real-space halo approach 

In this section, we briefly review the real-space dark matter halo approach (McClelland & Silk 1977; Scherrer & Bertschinger 1991; Sheth 
& Jain 1997; Yano & Gouda 1999; Ma & Fry 2000b,c), where the n-point correlation functions are described by the real-space convolution 
of halo density profiles. 

In the halo approach, we assume that all the matter is in halos with the density profile ph{x; m) parameterized by a mass rn. The halo 
mass is given by 



1 The nonlinear bispectrum has been well studied based on extended perturbation theory (Scoccimarro & Friemann 1999; also see Bernardeau et al. 2002a) 
or the recently developed dark matter halo approach (Ma & Fry 2000c; Scoccimarro et al. 2001; Cooray & Hu 2001a). 
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where VW is the virial volume of the halo. It is convenient to introduce the normalized halo profile defined as u m (x; m) — ph{x; m)/m, 
satisfying the condition 

d x Um(x; m) = 1. (2) 



In this paper we employ the virial volume as the boundary of a given halo, which can be formally defined, for example, based on the 
spherical top-hat collapse model. However, since in reality some matter is distributed outside the virial region, it is non-trivial to choose 
which boundary condition to use for the halo model. We will discuss how this uncertainty affects the halo model predictions, and show that 
it does not matter in the non-linear regime. 

The density field at some point can be given by the superposition of each halo density profile: p{x) = J^i ^"m, (a;; mi). The mean 
density of the universe is thus obtained from the ensemble average (Scherrer & Bertschinger 1991) 



po — (p(x)} — raju(xj — x; m,)| = Jdmmn(m) Jd 3 x'u m (x — x), (3) 

where we have replaced the ensemble average by a spatial average and an average over the halo mass function n(m): 
(Y]. S D (m — mi)8%(x' — a;,)) = n(m). Note that from equation (2) the integral over x' is unity, so the above equation is a statement 
about the normalization of the halo mass function. Throughout this paper we work in the comoving coordinate. 

Within the framework of the halo approach, the two-point correlation function (2PCF) of the density field can be expressed as the sum 
of correlations within a single halo (1-halo term) and between different halos (2-halo term): 

£(r) = ei h (r)+6a,(r) 

+ / dmi n(mi)— /rfm2ii(m2)— / dxiu mi (x — xi) / d X2 u m2 {x — X2)£,{xi — x% mi, m 2 ), (4) 
J Po J Po Jv vir Jv wit 

where £(xi — X2; mi,m<2) is the 2PCF of two halos with mass mi and 7712, and we have used r = x — x' . Here we have again employed the 
virial volume for the integration range. In the other words, the integration cutoff reflects the fact that we only account for mass contributions 
out to the virial region. If we require that £(r) is a function of the separation r alone (from statistical homogeneity and isotropy), then the 
density profile u m which we need here is an average over all halos of a given mass so that the resulting 2PCF has no particular direction 
dependence of r (see similar discussion in Seljak 2000). This holds for a spherically symmetric density profile, which is assumed in this 
paper for simplicity. If one considers non-spherically symmetric profiles, it requires some extra averages over all possible shapes of halos of 
a given mass. Equation (4) shows that we can obtain £ih(r) by a 3-dimensional integration for a spherically symmetric density profile, while 
the 2-halo term needs an 8-dimensional integration at most. 

In the same spirit as the derivation of equation (4), we can derive expressions for the 1-halo contribution to the n-point correlation 
function (n > 2): 



(5) 



Interestingly, this equation means that we can obtain any n-point correlation function by a 4-dimensional integration, once the density profile 
and the mass function are given. The 1-halo contribution is expected to be dominant in the strongly non-linear regime (<5 >• 1), and so the 
real-space halo approach can be useful for calculations in this regime. 

In this paper, we focus on the halo model predictions of the 3PCF, which is the lowest order statistic for describing non-Gaussianity of 
the density field. The 1-halo term can be expressed as 



(6) 



Under the assumption of statistical symmetry, the 3PCF can be described by three independent parameters characterizing triangle configu- 
rations. Likewise, the 2- and 3-halo terms in 3PCF can be derived within the real-space halo approach (Ma & Fry 2000c); for example, the 
2-halo term becomes 

C,2h(ri,r 2 ,r i ) = dmn(m)\™) dm' n(m')^- 
J \PoJ J Po 

f f 

x / d A x u m (x + ri)u m {x + r 2 ) / d^x' u m i {x + r 3 )£(\x — x'\ ; m, m) + perm(l, 2, 3). 

Jv vir Jv . 



(7) 



This equation implies that we have to perform at most an 8-dimensional integration to get ( 2 h for a given halo-halo correlation function, 
£(r; m, m'), while the 3-halo contribution to ( similarly requires the 12-dimensional integration. 
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2.2 Fourier-space halo approach 

Most recent studies of the halo approach has relied on the method build in the Fourier space (Seljak 2000; Ma & Fry 2000c; Peacock & Smith 
2000; Scoccimarro et al. 2001; Cooray & Hu 2001a,b; TJ02), since the Fourier-transformed counterparts of the n-point correlation function 
(power spectrum, bispectrum and etc.) can be expressed simply. 

Following the notation introduced in Cooray & Hu (2001a), the power spectrum can be expressed as the sum of the 1-halo and 2-halo 
contributions 

P(k) = P lh (k) + P 2h (k) = I? 2 {k,k) + [Il(k)] 2 P L (k), (8) 
with 



lfl(ki,...,kn) = J dm n(m)b r \m) (^j^ u m {ki) . . . i m (fcj. 



(9) 



Here u m (k) is the Fourier transform of the density profile and the following definition has been used recently (e.g., Seljak 2000; Scoccimarro 
et al. 2001; Cooray & Sheth 2002): 

l (k)= I 4irr 2 dru m (r)j {) {kr), (10) 



Jo 



where jo is the zeroth-order spherical Bessel function, jo(x) = sinx/x, and r v j r denotes the virial radius. In the following, quantities with 
tilde symbols denote the corresponding Fourier transformed quantities. Notice that all relevant quantities also depend on redshift, although 
we often omit it in the argument for simplicity. In contrast to the truncated halo profile above, Ma & Fry (2000c) employed the non-truncated 
profile whose Fourier transform is defined over an infinity integration range as u m (k) = J 4wr 2 drUm{r)jo(kr). The different halo 
boundary leads to different results in the Fourier-space halo model. In particular, we will focus on the difference between the predictions for 
the mass bispectrum in Ma & Fry (2000b), Scoccimarro et al. (2001) and Cooray & Sheth (2002) in comparison with the real-space halo 
model for the 3PCF in this paper. In the derivation of equation (8), we have assumed that the power spectrum for the halo-halo 2PCF can be 
given by P(k; mi, m 2 ) — b(mi)b(rn 2 )P L (k), where P L (k) is the linear power spectrum and b(m) denotes the bias parameter between 
the distribution of halos and the underlying density field. Through the definition £(r) = J k 2 dk/(2-K 2 )P(k)jo(kr), we can compute the 
2PCF of the density field from the power spectrum (8) within the framework of the Fourier-space halo model. The two derivations of £(r) 
based on the real-space and Fourier-space halo approaches should be equivalent to each other. However, the relation between £(r) and P(k) 
uses a Fourier transform over the infinite volume, while the integrations used in the halo model (see equations (4) and (10)) are confined to 
the virial volume of halos. This could lead to some discrepancies between the two halo model predictions for which will be carefully 
investigated below. 

Likewise, the Fourier-transformed counterpart of the mass 3PCF, the bispectrum, can be expressed as the sum of the three contributions 
where three dark matter particles reside in one halo (1-halo) or in two (2-halo) and three (3-halo) different halos, respectively: 

B(fci,fc 2 ,fc 3 ) = l!{ki,k 2 ,k 3 )+ [/ 2 1 (fci,fc 2 )J 1 1 (fc 3 )P i (fc3)+perm(l,2,3)] + J 1 1 (fc 1 )J 1 1 (fe 2 )/ 1 1 (fc 3 )B PT (fe 1 , k 2 , fe 3 ). (11) 

Here B PT denotes the perturbation theory bispectrum (e.g., see Jain & Bertschinger 1994) given by 



B PT (k u k 2 ,k 3 ) = 



10 , (k x | k 2 \ fci ■ fc 2 | 4 (k 1 ■ fc 2 ) 2 



7 \fc 2 kiJ k\k 2 7 k\k 



P L (k 1 )P L (k 2 ) +perm. (12) 



Here we have neglected the weak dependences of B on the cosmological parameter fi m0 . Using the bispectrum, the 3PCF can be defined 

as 

f 3 3 

C(ri,r 2 ,V3) = y IT (27r)' 3 ^ fcl ' fc2 ' fc3 ^ cxp [*( fel ■ r i + k 2 ■ r 2 + k 3 ■ r3)](27r) 3 <5i)(fci23), (13) 
where fci 23 = fci + k 2 + k 3 . 



2.3 Computational ease in real-space versus Fourier-space 

Before specifying the details of our method, we briefly compare the real-space and Fourier-space halo approaches for predicting the 3PCF of 
the mass distribution. 

Let us first consider the 1-halo contribution to the 3PCF. Equation (6) shows that, once the halo profile and mass function are given, the 
real-space halo model enables us to obtain the 1-halo term by a 4-dimensional integration for arbitrary triangle shapes. On the other hand, 
equations (11) and (13) imply that a 7-dimensional integration at most is required to get the 1-halo term based on the Fourier-space halo 
model, which is intractable with current computational expenses. Therefore, the real-space halo approach has a great advantage in predicting 
the 1-halo term of the 3PCF. 

However, the Fourier-space approach becomes useful for predicting the 2-halo and 3-halo terms, because the real-space model requires 
at most 8- and 12-dimensional integrations for calculating the 2- and 3-halo terms for given 2PCF and 3PCF of the halo distribution. On the 
other hand, from equations (11) and (13), one can see that the Fourier-space model allows us to obtain the 2- and 3-halo terms basically by 
a 7-dimensional integration. However, direct integration is still intractable. Therefore, in this paper we employ the Fourier-space model and 
develop approximations for calculating the 2- and 3-halo terms. Those approximations can significantly reduce the computational expenses 
and make the calculations tractable by regular numerical integrations, but still have adequate accuracy for our purpose. 
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2.4 Ingredients of the halo model 



To complete the halo model approach, we need suitable models for the three ingredients: the halo density profile, the mass function of halos 
and the biasing of the halo distribution, each of which depends on halo mass m and redshift z and have been well studied in the literature. 
We consider density profiles of the form 

Ps 



ph(r) 



(cr/r vil ) a (l + cr/r viT ) 3 



(14) 



where p 3 is the central density parameter and c is the concentration parameter. In most parts of this paper, we will take the NFW profile with 
a — 1 (Navarro, Frenk & White 1997). However, since other simulations with higher spatial resolution have indicated a w —1.5 (Fukushige 
& Makino 1997; Moore et al. 1998; Jing & Suto 2000), we will consider the effects of the variation in a on the 2- and 3-point correlation 
functions. The parameter p s can be eliminated by the definition (1) of halo mass as 



f r vir 



4nr 2 drp h (r) 



r 



2.Fi(3 — a, 3 — a, 4 — a, — c), 



(NFW; a = 1), 
(otherwise) , 



(15) 



where/ = l/[ln(l + c) — c/(l + c)] and2-F\ denotes the hypergeometric function. We employ the virial radius given by the spherical top-hat 
collapse model: m = (47rr^ ir /3)poA(z), where A(z) is the overdensity of collapse given as a function of redshift (e.g., see Nakamura & 
Suto 1997 and Henry 2000 for a useful fitting formula). We have A « 340 for the ACDM model. 

To give the halo profile in terms of rn and z, we further need to express the concentration parameter c in terms of m and z; however, 
this still remains somewhat uncertain. We consider the following c parameterized by the normalization and the mass slope: 



c(m, z) = co(l + z) 



m*(z = 0) 



(16) 



where m* (z — 0) is the nonlinear mass scale defined as S c (z = 0)/ D(z = 0)a(m*) = 1. Here S c is the threshold overdensity for spherical 
collapse model (see Nakamura & Suto 1997 and Henry 2000 for a useful fitting function), a(m) is the present-day rms fluctuations in the 
matter density top-hat smoothed over a scale R m = (3m/47rpo) 1,/3 , an d D(z) is the growth factor (e.g., see Peebles 1980). The values 
co ~ 10 and f3 ~ O(10 _1 ) are theoretically expected. The redshift dependence (1 + z) -1 is taken based on the numerical simulation results 
in Bullock et al. (2001). We take co = 10 and j3 = 0.2 as our reference model for the NFW profile, because TJ02 pointed out that the model 
can reproduce simulation results for the variance, skewness of weak lensing fields and the shear kurtosis. We will consider the influence of 
variations in the concentration parameter on the halo model predictions of the 2PCF and 3PCF. 

We will often use the following useful form for the Fourier transform of the NFW profile (Scoccimarro et al. 2001), which is needed for 
the Fourier-space halo model calculations: 



Um(k) = f 



sin v {Si(r?(l + c)) - Si(r?)} + cost? {Ci(r?(l + c)) - Ci(j?)} 



sin(?7c) 
»?(l + c) 



(17) 



where rj = fcrvir/c, Ci(x) = — dtcost/t is the cosine integral function and Si(x) = f*dtsm(i)/t the sine integral. The profile (14) 
with general a has the asymptotic behaviors of w m (fe) ~ 1 and oc k~ i+a for k <JC c/r v i r and k ^> c/r v i r , respectively. It is worth stressing 
that the equation above is derived from the NFW profile truncated at the virial radius to maintain mass conservation (see equations (1) and 
(3)). The truncated profile leads to the consequence that the 1-halo term of the power spectrum behaves like shot-noise as Pih(k) oc k° at 
small k. If one employs the non-truncated NFW profile, the resulting u m (k) is logarithmically divergent at k — > (see equation (22) in Ma 
& Fry 2000b) originating from the mass divergence from j°°4Tvr 2 drph(r), which does not yield the shot-noise power spectrum. 

For the halo mass function, we adopt an analytical fitting model proposed by Sheth & Tormen (1999), which is more accurate on cluster 
mass scales than the original Press-Schechter model (Press & Schechter 1974). The number density of halos with mass in the range between 
m and m + dm is given by 



n(m)dm 



m 

^A[l + (av)- p - 



auexp 



/ av\ dv 



where v is the peak height defined by 



5c(z) 



D(z)a(m) 



(18) 



(19) 



and the numerical coefficients a and p are empirically fitted from iV-body simulations as a = 0.707 and p = 0.3. The coefficient A is set by 
the normalization condition J°°duf(u) = 1, leading to A ^ 0.129. Note that the peak hight v is specified as a function of m at any redshift 
once the cosmological model is fixed. 

Mo & White (1996) developed a useful formula to describe the bias relation between the halo distribution and the underlying density 
field. This idea has been improved by several authors using TV-body numerical simulations (Mo, ling & White 1997; Jing 1998; Sheth & 
Lemson 1999; Sheth & Tormen 1999); we will use the fitting formula of Sheth & Tormen (1999) for consistency with the mass function (18): 



■ , x , av — 1 
b(u) = 1 + — ^ + 



2p 



5 C (1 + (au)p) 



(20) 
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where we have assumed scale-independent bias and neglected the higher order bias functions (62, 63, • • •)• This bias model is used for 
calculations of the 2-halo term in the 2PCF and the 2- and 3-halo terms in the 3PCF (see equations (8) and (11)). It should be noted that 
the requirement that the 2-halo term in the power spectrum gives the linear power spectrum in the limit k — > (w m (k) — > 1) imposes the 
condition jduf{u)b{u) = 1 for the integration mass range, since we have assumed that all the matter is in the form of virialized halos. 



3 ANGULAR iV-POINT CORRELATION FUNCTION 

In this section, we present the angular n-point correlation functions in analogy with the real-space halo approach to the 3D n-point correlation 
functions in §2. 1 . The results are applicable for a general projected cosmological field such as the angular galaxy distribution and weak lensing 
fields. 



3.1 Projected density field 

We first consider the 2-dimensional projection E(0) of the density field along the line of sight: 

E(0) = [d X W( X )p(x,d A 0), (21) 



where W(\) is me weight function, and \ and g(a(x) are tne comoving distance and comoving angular diameter distance, respectively. Note 
that x is related to redshift z via the relation d\ = dz/H(z). In this paper we take the projected field to be the weak lensing convergence 
field k (e.g., see Bartelmann & Schneider 2001) with the weight function given by 

W( X ,x.)--n m0 fl a , (22) 

where we have considered a single source redshift z s , corresponding to the comoving distance \s, for simplicity. The expressions derived 
below can be applied to the angular clustering of galaxies as well. Setting W(x) = dN/dx, the selection function of galaxies, and replacing 
the density field p by p ga iaxy (see §5 below for explicit expressions in terms of the halo occupation number of galaxies) will give the projected 
density of galaxies E ga i axy from equation (21), and can be used to replace E m in the equations below. 



3.2 The 1-halo contribution to angular n-point correlation functions 



Analogous to the 3D real-space halo approach, we can express the 1-halo contribution to the angular 2-point correlation function of the 
convergence field as 



w lh (6) 



-f 

Jo 



d (p £ m (<^; x)E m (y + 0; x), 



(23) 



where d 2 V/dxdQ is the comoving differential volume given by d 2 V/dxdil = d\(x) f° r a flat universe and we have assumed a circularly 
symmetric profile E m for halo of a given mass. We have also employed the flat-sky approximation, the Limber approximation (Kaiser 1992), 
which means that we account for contributions to w lh from lens structures at the same redshift, and the Born approximation. Finally, as in 
equation (4), we have restricted the integration range to the circular area Sl vir enclosed by the projected virial radius. Equation (23) implies 
that we can obtain win by a 4-dimensional integration if we have an analytical expression for E m . In the halo model, E m is the normalized 
column density field for a halo with mass rn defined as 



E m (6>) 



/ r vir 
dr\\u m {r\\, 



d A 0). 



(24) 



The factor d 2 V /dxdfl = d\(x) m equation (23) arises from the substitution d 2 x± = d 2 A d 2 <p, if one wishes to relate equation (23) to 
equation (4) for the 3-dimensional case. It is worth noting that, in analogy with equation (2), E m should be defined so that it satisfies the 
mass conservation condition inside the virial radius: 



/ d 2 r A 



£ m (0;x) = 4(x) / d 2 0E m (0) 



1. 



For the NFW profile (14), E m (0) can be analytically derived using formulae (2.266) and (2.269.2) in Gradshteyn & Ryzhik (2000) as 



(25) 
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where 6 v \ r is the angular size of the virial radius of a halo at a given distance X : $vir = r v i T /d A (x)- Note that our expression for S m differs 
from the result of Bartelmann (1996) because he took the projection over ry = [—00, 00] (his expression includes typos, so, e.g., see Wright 
& Brainerd 2000 for the correct expression). The difference turns out to be pronounced, e.g. using the expression for c — > 00 in the above 
equation leads to discrepancies between the real- and Fourier-space halo model calculations of w lh (9), as shown below. 

It is useful to replace the projected halo profile S m in equation (23) by the convergence field, n m (0), for a halo of mass m: 

K m {0; X ) = W( X ,Xs) E m (fl; X ) = ^Ga' 1 ~ X) f "'dry p h (r h d A ( X )0;m), (28) 

where we have used (3/2)f2 m o.ffo/po = 47rG. Hence, equation (23) can be rewritten in a more physically transparent form 

w lh {6)=J d X J dm n{m; x) J d 2 ip K m {(p;x)Km(<P + 0\x)- (29) 

This equation means that the 2PCF of the weak lensing convergence can be expressed as the line-of-sight integration of the lensing contri- 
bution of halos at different redshift, weighted with the halo number density. An interesting application of this formulation is the study of the 
2PCF of the reduced shear field, g = 7/(1 — k) (7 is the shear due to lensing) which is a direct observable of the cosmic shear measurement. 
Equation (29) no longer relies on the power spectrum to compute the 2PCF. Therefore, replacing n m in the above equation with g m for a 
given halo will yield predictions of the reduced shear 2PCF. So far cosmological interpretations of cosmic shear data have been made by 
comparing the data with the theoretical model of the shear 2PCF computed from a model of the 3D power spectrum (e.g., see Van Waerbeke 
et al. 2001). The non-linear correction in using the reduced shear rather than 7 could be important on sub-arcminute scales, where k is of 
O(10 _1 )for massive halos. This study will be presented elsewhere. 

Performing the 2-dimensional Fourier transform of equation (29) yields the 1-halo contribution to the angular power spectrum of the 
convergence field, C K {1): 



C K ,ih(l) = J dx ~^xi J ' dmn(m;x)\K m {l)\ 2 , (30) 

where we have assumed that the Fourier transform confined to the virial region is well approximated as the 2D Fourier transform over an 
infinite integration range and k m (l) is the 2D Fourier transform of n m defined as 

k m (l)= 2-K9d0 Km(0)J O (Ol). (31) 

Jo 

Here Jo (a;) is the zeroth-order Bessel function. The expression (30) coincides with the form in Cooray et al. (2000) (see also the original 
formulation in Cole & Kaiser 1988). 

To derive C K (l), the usual approach is to model the 3D power spectrum and then employ Limber's equation (e.g., see TJ02 and see 
references therein) to get 

Cih(1) = Jd X W 2 {x)d- A 2 {x)p[k = ^y;x) • (32) 
Using the expression (8) of P(k) based on the Fourier-space halo model leads to the another expression of wih(0): 

w lh (0) = Jd X W 2 { X )d- A 2 {x) jdmn{m- X )(j^j J l ^\u m (k)\ 2 J»(8l), (33) 

where k = l/d A ( X )- The difference between the real-space and Fourier-space approaches in equations (23) and (33) is only in the order of 
the projection. Therefore, the two approaches should be equivalent to each other (see also discussions in Cooray et al. 2000; Cooray & Hu 
2001a), if the weight function is a smooth function with respect to redshift as required for Limber's approximation. We will check this by 
comparing the predictions from the two methods 

In a similar manner, we can derive the 1-halo term for the angular n-point correlation function based on the real-space halo approach as 

w? h (0 u 2 ,...,0 n ) = Jd X W n ( X ) Jdmn(m; X ) J d 2 <p E m (^ + i; X )^m(<P + 9 2 ; x) • • • S m ( V + 0„; X ). (34) 

This equation implies that we can obtain the 1-halo term of any n-point angular correlation functions by a 4-dimensional integration. This 
expression will be useful in investigating statistical properties of the projected field at non-linear angular scales. In this paper, we focus on 
the 3-point correlation function of the weak lensing convergence, given by 

Z lh (9 1 ,e 2 ,e 3 ) = d X d A ( x )W : \x) Jdmn(m-,x) (j^J J d 2 <p Y, m { V + 0v, X )^m{v + Or, x)^m(v + 3 ; x) 

= / dxd 2 A (x) dmn(m; X ) dtp K m (<p + 0i; x)K m (<p + 6 2 )n m (v + 3 ; X ). (35) 
Jo J Jn vit 

From statistical symmetry, Z can be given as a function of three independent parameters characterizing the triangle shape on the sky. 
On the other hand, in the conventional approach to the 3PCF, Z is expressed in terms of the 3-dimensional bispectrum as 
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Figure 1. The 2-point correlation function of the mass density field, f (r), at z = for ACDM. The solid curve is the prediction from the Peacock & Dodds 
fitting formula. The thick and thin dashed curves show the 1-halo contribution to f (r), calculated using the real-space halo approach with different boundary 
conditions. The thick and thin dotted curves are the 1-halo and 2-halo terms from the Fourier-space halo approach. The dot-dashed curve is the total prediction 
of the 1-halo plus 2-halo terms, where we have used the real-space (thick dashed curve) and Fourier-space (thin dotted curve) models for the 1-halo and 2-halo 
predictions, respectively. 



z(e 1 ,e 2 ,e 3 ) = jf S d x w 3 ( x )d^ 4 J f[ k 2 , k 3 - x ) ■ o^fsl^), 06) 

where k% = k/d>A(x)- We use this formulation for evaluating the 2-halo and 3-halo contributions to Z, using the expressions (11) for the 2- 
and 3-halo bispectrum. The resulting 6-dimensional integrations are intractable, so we develop approximations, presented in Appendix A, to 
evaluate them. 



4 RESULTS 

In this section, we present halo model predictions for the three-point correlation function for the mass and galaxy distribution and for the 
weak lensing field. We employ a scale invariant spectrum of primordial fluctuations with the BBKS transfer function (Bardeen et al. 1986). 



4.1 2-point correlation function 

Before considering the 3PCF, we demonstrate the validity of the real-space halo model by comparing our results for the two-point correlation 
function with those of the Fourier-space halo model well studied in the literature. 

For the NFW profile, the 1-halo term in equation (4) can be further analytically simplified as 

M r ) = ^Jdmn(m) (j^J ^£^j J^ds (1 + cs/r vir )- 2 l(s;r), (37) 
with 



0, (r > 2r v ir, or r — s\ > r v jr) 

l(s;r) = { r(l + c{r T -s\/r vil )(l+c)> (s + r > r vil ) (38) 

(y* _|_ g\ y g\ 

—7- 1 ^-r-. — —z-j] r-; r- , (otherwise), 

r(l+c\r-s\/r vil )(l + c(r + s)/r viI ) v ' 

where we have taken the virial volume for the integration range of d 3 s in equation (4) and the boundary condition u m (x) = for x > r v - a . 
The predictions for £1^ on nonlinear scales are not sensitive to the boundary condition, as shown below. Note that, although the integration 
over ds can be further analytically calculated as done in Ma & Fry (2000c) and in Sheth et al. (2001), the resulting expression for is 
lengthy for the boundary condition and so we stopped at the above expression. 

Figure 1 plots the mass two-point correlation function £(r) at z = for the ACDM model. The halo model prediction is shown by 
the dot-dashed curve, which matches the solid curve given by the fitting formula of Peacock & Dodds (1996; hereafter PD). We consider 
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Figure 2. (a) The left panel shows the dependence of the 1-halo contribution to f (r) on the maximum mass cutoff used in the calculation. From top to bottom, 
the maximum mass is 10 16 , 10 15 , 10 14 , 10 13 and 10 12 M Q . For comparison, the dashed curve denotes the full £(r) (lh+2h) for our reference model, (b) In 
the middle panel, the three solid curves demonstrate the dependence on the inner slope parameter a of the halo profile (see equation (14)). (c) The right panel 
shows the dependence on the concentration parameter. The result for our fiducial model (cq, /3) = (10, 0.2) is shown by the thick solid curve. 



the NFW profile with the concentration specified by (co, 13) = (10, 0.2) in equation (16) as the fiducial model. During the preparation of 
this paper, Smith et al. (2002) have proposed a new fitting formula of the non-linear power spectrum that can better match high-resolution 
iV-body simulations. For the ACDM model, the effect is small: the power spectrum differs from the PD prediction by ~ 10% at most at 
k ;$ 30/i/Mpc. The thick dashed curve shows the real-space halo model prediction for the 1-halo contribution to £ computed from equation 
(37), while the thick dotted curve is the prediction from the Fourier-space halo model. The two halo model predictions are in remarkable 
agreement for small scales r ^2 h~ 1 Mpc, but deviate slightly at the transition scale of ~ 2 h~ 1 Mpc between the non-linear and linear 
regimes. The discrepancy between the real-space and Fourier-space halo models could be ascribed to the boundary condition used in the 
integration. The thin dashed curve shows the result obtained without using the boundary condition (u m (x) = for x > r v i r ) in equation 
(37), implying that the change alters the predictions only on scales r ^ 1 ft~ 1 Mpc. However on these scales the 2-halo term dominates 
as shown by the thin dotted curve. Hence, to summarize the results in this figure, the real-space halo model can be used to predict the 
statistical properties of the non-linear density field with the same accuracy as the Fourier-space halo model whose validity has been carefully 
investigated in the literature (Seljak 2000; Ma & Fry 2000c; Scoccimarro et al. 2001; TJ02). 

The left panel of Figure 2 shows the dependence on the maximum mass cutoff used in the calculation on £(r). From top to bottom, the 
solid curves show the results for the maximum cutoff of 10 16 , 10 15 , 10 14 , 10 13 and 10 12 M Q . Massive halos with m > 1O 14 M yield the 
dominant contribution (greater than ~ 80%) to for r £ 0.1 h~ 1 Mpc, while less massive halos become relevant at smaller scales. The 

dashed curve denotes the full halo model prediction for £(r) (lh+2h) as shown in Figure 1. 

The middle and right panels of Figure 2 show the dependences of £(r) on the variations in the halo profile. From the middle panel, one 
can see that increasing a for the halo profile of equation (14) leads to a higher amplitude and steeper slope for £(r) on small scales (see also 
Jain & Sheth 1997; Ma & Fry 2000c). For generic a, the angular integration in J d 3 s in equation (4) is not analytic, so we performed the 
3-dimensional integration to get £ih (r) numerically. The right panel shows the dependences on the concentration parameter parameterized 
in terms of the normalization co and the slope (3 as c(m) = co(m/m„) - ' 3 . The result of our fiducial model (co, (3) = (10, 0.2) is shown by 
the thick solid curve. It is apparent that increasing Co with fixed f3 or decreasing /3 with fixed co increases £ at small scales. This is because 
these variations lead to more concentrated density profiles for halos more massive than the nonlinear mass scale m* . Since the massive halos 
dominate the contribution to £ at the scales considered here, this has the effect of increasing £ , The results in Figure 2 imply that we can 
adjust the inner slope a and the concentration parameter, both of which are somewhat uncertain theoretically and observationally, so that the 
halo model can reproduce the PD result for a given mass function (see Seljak 2000 for a similar discussion). We will argue that the 3PCF can 
be used in combination with the 2PCF to constrain properties of the halo profile. 

Figure 3 shows the 2PCF of the weak lensing convergence field. The real-space halo model prediction (dashed curve) for the 1-halo 
contribution is compared with the result of the Fourier-space halo model (thick dotted curve), computed from equations (29) and (33), 
respectively. Here we have considered the NFW profile and source redshift z s = 1 for simplicity. The different halo model predictions 
perfectly agree over the angular scales we have considered. This agreement is encouraging, because the real-space halo approach can then be 
used to predict the 3PCF with the same accuracy expected as for the Fourier-space model well studied in the literature (TJ02 and references 
therein). For comparison, the solid curve denotes the PD result and the thin dotted curve is the Fourier-space halo model prediction of the 
2-halo term. It clarifies that the 2-halo term is important on 9 ^ 3' for the ACDM model. The dot-dashed curve is the total contribution to the 
2PCF from the 1- and 2-halo terms, which slightly deviates from the PD result on the small scales 9 £ 1'. As discussed below equation (26), 
if we use the expression in Bartelmann (1996) for the NFW convergence field in the calculation of the 1-halo term, one finds discrepancies 
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Figure 3. Comparison between model predictions of the two-point correlation function of the weak lensing convergence field, as in Figure 1 . The dashed curve 
is the real-space halo model prediction (equation 29) for the 1-halo contribution to the 2PCF, while the thick dotted curve denotes the result of the Fourier-space 
halo model computed from equation (33). The top, broken curve shows the overestimation of the halo model result (lh+2h), if the NFW convergence field in 
our equation (29) is replaced by Bartelmann's (1996) result. 




Figure 4. The sketch of the triangle configuration parameterized by r, q and ip used in this paper to describe the 3PCF. 



between the real-space and Fourier-space halo models. The broken curve shows the result (1-halo + 2-halo terms), which overestimates w K 
by 10-25 % over the scales considered. 



4.2 3-point correlation function 

4.2.1 Dark matter correlation function 

We first consider the 3PCF of the mass density field. To obtain the 1-halo contribution, where all the three matter particles reside in one halo, 
we perform the following 4-dimensional numerical integration; 

I 1 / \ ^ /* r vir 1*277 /*7T 

Cih(r, q, ip) = Jdmn(m)i — \ J ds J dip J d9 s 2 sin 9 u m (s)u m (\s + r\)u m (\s + q\), (39) 

where |s + r\ = (s 2 + r 2 + 2sr cos 9) 1 ^ 2 and |s + q\ = [s 2 + q 2 + 2sg(sin ip sin 9 cos ip + cost/; cos O)] 1 ^ 2 . Note that we employ the 
boundary condition u m (r) — for r > r v - 1Y , but this does not affect the final results on the non-linear scales, as explained in Figure 1 and 
explicitly demonstrated in Figure 8. In the above equation, we have used the fact that the 3PCF can be expressed as a function of the three 
independent parameters r, q and ip specifying the triangle configuration (see Figure 4). 

To complete our halo model predictions, we develop approximations for calculating the 2- and 3-halo terms in Appendix A. The 2-halo 
term is relevant only for small range of scales in the transition between the quasi-linear and non-linear regimes (r ~ 1 /t -1 Mpc), while the 
3-halo term is dominant on larger scales. As in the literature, we consider the hierarchical 3PCF amplitude defined as 
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Figure 5. The halo model predictions for the reduced 3PCF, Q, as a function of the shapes of triangle configurations. The results for the NFW profile at z = 
are plotted for the ACDM model. The triangle shape is parameterized by three parameters r, q and ip as illustrated in Figure 4. The four panels show the 
results for Q vs. ip, with r = 0.02, 0.1, 1.0 and 10 h~ 1 Mpc as indicated in the panels. In each panel, the solid, dashed, dot-dashed, dotted, and broken curves 
denote the Q parameter for r/q = 1, 2, 3, 4 and 5, respectively. 



Q(r,q,i>) = - Ji r,q, 1L, ( iv (40) 

Z{r)Z{q) + Z(r)Z(\r - q\) + ?(g)£(|r - q\) 

where C, = Cih + ^2h + ^3h- In the following halo model predictions for Q, we will maintain self-consistency by using the halo model for 
calculating the 2PCF in the denominator of Q. 

Figure 5 plots the dependences of Q(r, q, ip) on the size and shape of triangle configurations. In this figure we have used the NFW 
profile. The four panels show the results for r — 0.02, 0.1, 1 and 10 /i _1 Mpc. In each panel, the solid, dashed, dot-dashed, dotted and 
broken curves plot the Q values for q/r = 1, 2, 3, 4 and 5 as a function of the interior angle tp. It is apparent that on the highly non-linear 
scales r £ 0.1 ft _1 Mpc the 3PCF has a weak dependence on the triangle configuration; even so, for plausible halo model parameters the 
hierarchical form Q w const, does not hold in that Q depends both on the size and shape of the triangle configuration (also see discussion 
in Ma & Fry 2000b). The curves with q/r = 1 show that Q decreases with decreasing tp at tj) ^ 0.2n because of the strong suppression due 
to £ 2 in the denominator of Q. On larger scales r ^ 1 Mpc, an oscillatory feature in Q appears as predicted by perturbation theory (e.g., 
Fry 1984b; Barriga & Gaztanaga 2002). Thus, the configuration dependence of the 3PCF is more prominent on quasi-linear scales than on 
strongly non-linear scales. 

The features observed in Q can be understood using Figure 6, which shows the 1-, 2- and 3-halo contributions to Q separately for the 
triangle configurations in Figure 5. Note that these separate contributions are shown for the 3PCF, the numerator of Q, while the 2PCF factors 
in the denominator include the full contribution from the 1- plus 2-halo terms. This figure clarifies that the 1-halo term yields the dominant 
contribution to the 3PCF on small scales r 1 Mpc. The 2-halo term becomes relevant over the transition scales 1 ^ r ^ 5 Mpc between 
the non-linear and linear regimes, while the 3-halo term dominates for r ^ 10 Mpc. 

Figure 7 displays the halo model predictions for Q for equilateral triangles as a function of the side length r. The three dashed curves 
plot the 1-, 2- and 3-halo contributions separately. For comparison, the perturbation theory result is shown by the dotted curve. The figure 
again clarifies that, at the small scales r ;$ 0.5 /i _1 Mpc, the 1-halo term dominates the total contribution, while the 3-halo term captures 
the large-scale correlations in the quasi-linear regime. A comparison between the perturbation theory result and the 3-halo term reveals that 
our approximation (A4) reproduces the perturbation theory result at large scales r ^ 5 /i _1 Mpc. Thus the way the halo model provides a 
separate description of the 1-, 2- and 3-halo terms can clarify how gravitational clustering transits from the linear regime to the strongly 
non-linear regime as one goes to smaller scales. 

The halo model predictions match TV-body simulation results on small scales r £ 1 ft -1 Mpc (R. Scoccimarro; private communication). 
However, the Q parameter from the halo model has a bump feature at r ~ 2 /i _1 Mpc which does not seem to be present in simulation data. 
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Figure 6. The 1-, 2- and 3-halo contributions to the Q parameter are shown by the dashed, dot-dashed and dotted curves, respectively. The triangle shapes are 
parameterized as in Figure 5. Note that in the top panel the 1-halo term dominates Q, while in the bottom panel the 3-halo term dominates. 



The corresponding bump feature in the bispectrum has been found at k ~ 1 /iMpc -1 in previous work (see Figure 2 in Scoccimarro et al. 
2001 and Figure 13a in Cooray & Sheth 2002). Hence, it is unlikely to be due to inaccuracies in our approximations used for the 2- and 
3-halo term calculations, although the approximations tend to overestimate the true amplitudes to some extent as carefully investigated in 
TJ02. It probably does not reflect real properties of dark matter clustering either (as discussed below for Figure 10). Rather, it can be ascribed 
to an inaccuracy in the standard implementation of the halo model which appears only on transition scales of order a few Mpc, in contrast to 
the success of the halo model in the non-linear and linear regimes. 

One possible origin of the inaccuracy is the sharp cutoff at the virial radius for the integrals for the 2PCF and 3PCF, as have been used 
for the real-space halo model in this paper as well as for the Fourier-space model in Scoccimarro et al (2001) and Cooray & Sheth (2002). As 
shown in Figure 1, modifications of the integration range could alter the halo model prediction at the transition scales. In the upper panel of 
Figure 8, the thick dashed curve shows the result obtained without using the boundary condition u m (x) = for x > r v j r for the calculations 
of the 1-halo terms of the 2PCF and 3PCF, which enter the denominator and numerator of Q, respectively. One can see that the bump feature 
is weakened through complex dependences - the 2PCF amplitude is enhanced by the modified boundary condition more strongly than the 
3PCF. In fact, this trend is verified by the results in Figure 3 and 4 in Ma & Fry (2000b), which show that the Fourier-space halo model 
yields bispectra with no bump feature. Ma & Fry employed the Fourier transform of the non-truncated profile, u m (k), calculated over an 
infinite integration range, which includes mass contribution outside the virial region. However, in this paper we have used the virial boundary 
condition since it preserves mass conservation. We have also clarified in Figure 3 that, if we include mass contributions beyond the virial 
radius, the halo model prediction for the lensing 2PCF overestimates the amplitude (see the dot-dashed and broken curves in Figure 3). In 
reality, the mass distribution outside the virial region eventually merges into quasi-linear structures, such as filamentary structures in the 
universe, which are relevant for the transition scale (S ^ 1), and are unlikely to follow the halo profile far outside the virial region. Therefore, 
a careful investigation of which boundary conditions to use for the halo model integrals will be needed to achieve more accurate predictions 
at the transition scale. 

Another effect that could be manifested at these scales is the exclusion effect of different halos for the 2- and 3-halo terms. We could 
impose the condition that different halos are separated by scales larger than the sum of their virial radii. The halo model does not explicitly 
account for this effect. Since the transition scale ~ 1 Mpc is of the order of the virial radius of massive halos that make a significant 
contribution to the integrals, the halo model used here might overestimate the contribution from the 2- and 3-halo terms at these scales. 
The dot-dashed curve in the lower panel of Figure 8 shows an estimate of this effect. The estimate is obtained by imposing the condition 
r v i T (m) < r/2, an approximate prescription to exclude pairs and triplets of halos contributing to the 2- and 3-halo terms at scales smaller 
than the sum of their virial radii. The plot shows that including the exclusion effect suppresses the bump feature in Q. The dotted curve in 
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Figure 8. This figure shows the effect of modifying the halo model calculation by altering the halo boundary condition and by exclusion effects in the 2- 
and 3-halo terms. The upper panel shows how the bump in Q C q(r) is suppressed if halos are not truncated at the virial radius (see text for details). The thick 
dashed curve shows Qoq(f ) and the three thin dashed curves show the 1-, 2- and 3-halo contributions. For comparison, the solid curves are the halo model 
results shown in Figure 7, in which a sharp cutoff is used at the virial radius. The dot-dashed curve in the lower panel shows the effect on Q C q(V) if halo 
exclusion effects are included in calculations. The 2-halo term of the 2PCF and the 2- and 3-halo terms of the 3PCF, which enter denominator and numerator 
in Qoq(i"), are constrained by the condition r v i r (m) < r/2. This effect suppresses the amplitudes of the 2- and 3-halo contributions to Q C q, again leading 
to the suppression of the bump feature. The 1-halo term is not shown in this panel as it remains unaffected. The thick dotted curve shows the prediction for 
Qoq(r) if both boundary and exclusion effects are included. 
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Figure 9. The results for Q C q(r) are shown for different cosmological models as in Figure 7. The dotted curve denotes the ACDM model with erg = 1.1, 
compared to ag = 0.9 for our fiducial model. The dashed curve shows the result for the SCDM model. To show the sensitivity of Q to the shape of the matter 
power spectrum, we also show the result for the rCDM model, which is different only in the shape parameter (r = 0.21) from the SCDM model (r = 0.5). 



the lower panel of Figure 8 shows the result of applying both alterations (boundary and exclusion effects) to the calculations. A systematic 
resolution of these problems is beyond the scope of our paper, and will be considered elsewhere. Finally, the figure implies that, except for 
the problematic bump feature, a transition between the non-linear amplitude (Q ^ 3) and the quasi-linear amplitude (Q 1) occurs at scales 
of a few /i _1 Mpc, corresponding to a flattening of the 2PCF £(r) (see Figure 1). 

Figure 9 shows the dependences of Q on the cosmological model. We consider three models that differ from our fiducial model. One is 
the A CDM model normalized with as = 1.1, motivated by the possible detection of the Sunyaev-Zel'dovich effect in the CMB (Bond et al. 
2002). The second is a flat model without cosmological constant (SCDM), with parameters f2 m o = l,h = 0.5 and as = 0.6. To clarify the 
dependence on the shape of the matter power spectrum, we also consider the rCDM model, where the shape parameter T — 0.21, but the 
other parameters are same as in the SCDM model. The figure reveals that the Q parameter is not very different for the different models but 
does vary slightly with changes in the various model parameters. Thus accurate measurements must be interpreted with precise predictions 
that explore all parameter dependences. 

For completeness we present our results with the triangles specified using a set of alternative parameters which have been used in the 
literature (e.g., Peebles 1980): 



T23 
ri2 ' 



r-ji — r 2 3 
ri2 



(41) 



with the condition ri2 < r^z < rz\ which imposes the constraints u > 1 and < v < 1. The 3PCF for any triangle can then be expressed 
as a function of three parameters ri2, it and v, where it and v characterize the shape and ri2 the size of a triangle. 

Figure 10 shows the dependence of Q on ri2, u and v. The values of r 12 , u and v are chosen such that our model predictions can 
be compared with the TV-body simulation results shown in Figure 2 in Jing & Bomer (1998). The data with the error bars represent the 
simulation results, which were kindly made available to us by Yipeng Jing. Although the cosmological parameters in Jing & Borner (1998) 
are slightly different from our model, this fact does not matter because of the weak dependence of Q on the cosmological model as explained 
in Figure 9. It is also noted that the halo model prediction does not account for the effect of averaging over bin widths in it and r\i which 
is taken for the Q measurement from the simulations. This figure reveals that our halo model agrees well with the simulation results for 
the configuration dependence of Q and its amplitude at the non-linear scales ri2 = 0.2 and 0.32 /i _1 Mpc, as well as the quasi-linear scale 
ri2 = 3.25 /i _1 Mpc. However, as discussed above, for the case in which one side length of the triangle is comparable to 1 /i _1 Mpc, the halo 
model tends to overestimate the simulation result. More detailed comparisons with recent higher resolution simulations will be presented 
elsewhere. 

Figure 11 shows how the 1-halo contribution to the mass 3PCF depends on the maximum mass cutoff used in the calculation against 
the side length of equilateral triangles. This can be compared to the dependence of the 2PCF on mass cutoff shown in the left panel of Figure 
2. For the 3PCF more massive halos dominate the contribution at a given length scale: for example, at r = 0.5 h~ Mpc, over half the 
contribution to the 3PCF is from halos with m > 1O 15 M0, whereas for the 2PCF the mass range is m > 1O 14 M0. 

In Figure 12 and 13, we present dependences of the Q parameter on possible variations in the halo inner slope and the concentration 
parameter, as in the middle and right panels of Figure 2. To do this, we use the parameters r, q and if; for the triangle configuration as in 
Figure 5. Note that the results shown are computed from only the 1-halo contributions to the 3PCF and 2PCF, which enter the numerator and 
denominator of Q, respectively, since the variations of the halo profile are relevant only at the non-linear scales. Increasing the inner slope 
parameter a for the profile (14) leads to higher Q (see similar discussions in Ma & Fry 2000a for the bispectrum). A comparison of Figures 
12 and 13 shows that Q is more sensitive to a for triangles with smaller size. Compared with the effect on the 2PCF shown in Figure 2, this 
result also means that the 3PCF is more sensitive to the inner slope of halo profile than the 2PCF. We find the following fitting formula of the 
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Figure 10. Halo model predictions for the dependence of the Q parameter on the triangle shapes parameterized by ri2, u and v (see equation (41)). The values 
of ri2, u and v are the same as those used in Figure 2 in Jing & Borner (1998), where the Q parameters were computed from ./V-body simulations. The data 
with error bars denote the simulation results. 




r (h 1 Mpc) 



Figure 11. The dependence of the mass 3PCF on the maximum mass cutoff used in the halo model calculations, plotted against the side length of the equilateral 
triangle. From top to bottom, the maximum mass is 10 16 , 10 15 , 10 14 , 10 13 and 10 12 Mq. The dashed curve is the full 3PCF. 
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Figure 12. The dependences of the Q parameter on the halo profile parameters: the inner slope a and the concentration parameter, c(m), parameterized 
in terms of the normalization co at the nonlinear mass scale, and the slope /3. The upper panels show the results for Q versus if) with the two side lengths 
r = q = 0.02 h _1 Mpc (Figure 4 defines ip, r, q). The lower panels have r = q/3 = 0.02 h~ Mpc. 



dependence of the 2PCF and the 3PCF on a, applicable over the range 0.5 < a < 1.5; 

,2/ n r ^2 1 l+0.017(a-l)+0.074(a 2 -l) \ > l+0.021(a-l)+0.0S(a 2 - 1) „ ^ l+0.094(a-l)+0.24(a 2 -1) 

? \ r ) = l^NFwl , ((r)=(,NFW Q = Qnfw (42) 

for r = 0.02 /i _1 Mpc, while 

fir) = [^F W ] 1+0 - 048(Q - 1)+0 -° 34(q2 - 1) , CM = Cjvfw 1+0-^("-1)+0.033( Q 2 -1 )i q = Qjvfw 1 + 0.38( Q -1) + 0.057( q 2 -1) ^ 

for r = 0.1 /i _1 Mpc. Here £nfw, Cnfw and Qnfw denote the use of the NFW profile (a = 1); we have used equilateral triangle 
configurations for the 3PCF calculation. Figure 14 explicitly plots the a dependence of Q. 

The right panels in Figure 12 and 13 illustrate the sensitivity of Q to possible variations in the halo concentration parameter (16). One 
can see that the change in concentration affects the Q amplitudes, but does not strongly alter the configuration dependence. Increasing Co 
with fixed f3 or decreasing f3 with fixed Co leads to a higher amplitude for Q, which is the same trend as for the 2PCF shown in Figure 2. This 
implies that the 3PCF is more sensitive to the halo concentration than the 2PCF. These results in Figure 12 and 13 show that measuring the 
3PCF could constrain halo profile properties complementary to the measurements of the 2PCF, in analogy with determinations of the galaxy 
bias parameter from measurements at quasi-linear scales (Feldmann et al. 2001; Verde et al. 2002). 

In Figure 12-14, we have demonstrated the sensitivity of the mass 3PCF to the halo profile properties. We have compared our halo model 
predictions for Q with the asymptotic shape derived in Ma & Fry (2000b) based on the halo model. The asymptotic Q amplitude depends 
only on the mass slope /3 of the concentration parameter (see equation (16)), the slope a' in the low mass end of the mass dependence of the 
mass function, n(m) oc v a , and a primordial spectral index n. From equation (7) in Ma & Fry (2000b), the asymptotic Q for equilateral 
triangles is given by 

^/ \ -i a'(n + 3) 

Q{r)(xr > ^2WTTy (44) 

where /3' w 0.8/3 and we have considered equilateral triangle configuration. This equation was derived under three assumptions: the con- 
tribution from the exponential cutoff of the mass function in the high mass end was ignored, a scale-free primordial power spectrum was 
employed and the fc-dependence of the Fourier-transformed halo profile, u m (k) was ignored. The third assumption leads to the consequence 
that the asymptotic Q has no dependence on the inner slope of the halo profile. Since for the profile in equation (14) u m (k) does depend on 
the inner slope as u m {k) w 1 and oc k~ +a for k <§C c/r v i r and k 3> c/r v i r , the asymptotic Q should depend on the inner slope or more 
generally the halo profile shape. By comparing our predictions of Q with the analytical result of Ma & Fry (2000b), we find that our halo 
model predicts a steeper slope for Q(r) than expected from the asymptotic shape in equation (44), for plausible values of j3 and a' . 

This discrepancy is probably due to the assumptions used in analytically deriving the asymptotic shape. The r-dependence of Q arises 
from a complex superposition of contributions from the halo profile, the power spectrum shape and the shape of the mass function. Indeed, 
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Figure 13. The Q values, as in the previous figure, for r = 0.1 h 1 Mpc. 
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Figure 14. The Q parameter for equilateral triangles as a function of the inner slope parameter a for the halo density profile given by equation (14). The solid 
and dashed curves are the results for r = 0.02 and 0.1 fe _1 Mpc, respectively. 



as an example of this possibility, Taruya, Hamana & Kayo (2002) explicitly showed, based on the halo model, that the non-Gaussian tail of 
the probability distribution function of the density field arises from such a superposition for small smoothing scales. 

An important implication pointed out in Ma & Fry (2000b,c) is that plausible halo model properties are unlikely to follow the stable 
clustering hypothesis. This hypothesis is useful because it allows us to analytically predict the behavior of non-linear gravitational clustering 
(e.g. Peebles 1980; Jain 1997). In fact, the popular PD fitting formula widely used in the literature was derived based on this hypothesis (also 
see Hamilton et al. 1991; Jain, Mo & White 1995). Very recently, Smith et al. (2002) showed that the non-linear power spectra measured 
from the high-resolution TV-body simulations indeed showed a weak violation of the stable clustering hypothesis. Figure 9 in Smith et al. 
(2002) compares the fc-slope of the measured non-linear power spectrum with the asymptotic prediction from Ma & Fry (2000b), which also 
shows some discrepancies. Hence, a careful investigation will be needed to clarify how non-linear gravitational clustering can be described 
by the halo model ingredients in connection with the stable clustering hypothesis and the hierarchical ansatz (see also the next subsection). 
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Figure 15. The hierarchical amplitude Q n for the mass n-point correlation function in the strongly non-linear regime is shown versus the order n. We consider 
the configuration with n equal sides of length r = 0.1/i~ 1 Mpc and equal interior angles. The solid curve shows the halo model prediction. The dashed curve 
shows the scaling for Q n as a function of Q3 and n proposed by Fry (1984a), while the dotted curve shows the formula proposed in Hamilton (1988) (we 
have used the halo model result for Q3 in these formulae). The dot-dashed curve is the prediction from hyper-extended perturbation theory (Scoccimarro & 
Frieman 1999). 



4.2.2 The hierarchy of higher-order correlation functions 

As stated in §2.1, the real-space halo model can be a useful tool for predicting the n-point correlation functions in the strongly non-linear 
regime, where correlations within one halo dominate the contribution to the n-point functions. The exact expressions for the the hierarchy 
are provided by the BBGKY equations (Peebles 1980), which govern collisionless gravitational system. However, their complex form makes 
it intractable to obtain solutions for the higher order correlations in the non-linear regime. Here we will use the halo model to shed some light 
on a fundamental question about gravity that has been discussed in the literature over the last decades: what is the asymptotic small-scale 
behavior of the hierarchy of the n-point functions under gravitational clustering (e.g., Peebles 1980; Fry 1984a,b; Hamilton 1988; Ma & Fry 
2000b; Bernardeau et al. 2002a)? 

The reduced n-point correlation function, Q n , is defined as 



Qn 



(«), 



ri, 







SlabcliiigsH edgea ij^(rij) 



(45) 



where the denominator is given by all topologically distinct tree diagrams (the different n™~ 2 ways of drawing (n — 1) links that connect n 
points). This form is motivated by the expectation oc £ n_1 , indicated observationally by the pioneering measurements of the two- and 
three-point functions from galaxy surveys (e.g., Groth & Peebles 1977; Peebles 1980), and theoretically by perturbation theory (Fry 1984b). 
The amplitudes Q n are a natural set of statistics to describe the non-Gaussianity that results from gravitational clustering. 

The solid curve in Figure 15 shows the halo model predictions for the hierarchical amplitudes Q n for 3 < n < 6. We consider 
configurations with n equal sides, each of length r = 0.1/i _1 Mpc, and equal interior angles for evaluations of the n-point functions. We 
take the halo profile to have the NFW form. For comparison, the dashed curve shows the scaling law for Q n proposed in Fry (1984a): 
Q n = (4Q 3 /n) ( "" 2) n/(2n - 2). The dotted curve shows the scaling suggested by Hamilton (1988): Q n = (Q 3 /n) n " 2 n!/2. The dot- 
dashed curve is the prediction from hyper-extended perturbation theory for 3 < n < 5 (Scoccimarro & Frieman 1999). The halo model 
predicts stronger clustering with increasing n than the formulae of Fry and of Hamilton, but not as strong as hyper-extended perturbation 
theory. Note that higher-order correlations are more sensitive to variations in the inner slope and halo concentrations of massive halos. 

The results in Figure 15 might be compared with Figure 3 in Szapudi et al. (1999), which show the cumulants S n (3 < n < 10) 
measured from the iV-body simulations. From equation (45) we can expect following rough relation between the cumulants and the reduced 
n-point correlation functions: 



S n (r) 



K (a)( r )](»-i) 



Q n (r). 



(46) 



Table 1 shows a comparison of the simulation results for S n in Szapudi et al. (1999) with the theoretical predictions scaled from the values 
of Q n in Figure 15. Note that the simulation results were read from the dashed curves in Figure 3 in Szapudi et al. (1999) for the SCDM 
model (the cosmological parameters are slightly different from the model we have considered), and we expect similar values for the ACDM 
because of the weak dependence of S n on cosmological models. One can see that the rapid increase in the amplitudes of S n with n seen in 
the simulation result is better described by the halo model predictions than those of Fry and of Hamilton. However, as explicitly pointed out 
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Table 1. Comparison of the theoretical predictions for the cumulants, S n , at r = 0.1/i _1 Mpc with the simulation result. The analytical predictions are 
estimated from the reduced rt-point functions, Q n , shown in Figure 15 using equation (46). The simulation results correspond to the values denoted by the 
dashed curves in Figure 3 in Szapudi et al (1999), which are measured from A^-body simulations for the SCDM model. 
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in Szapudi et al. (1999), the simulation results might be overestimated due to the shot noise effect. Therefore, further detailed work is needed 
to systematically explore the consistency of the halo model with gravitational dynamics. 



4.2.3 The 3 -point correlation function of galaxies 

It is straightforward to extend the dark matter halo approach to calculate the 3PCF of galaxies which can be measured from galaxy surveys. 
To do this, we need to know the mean number of galaxies per halo of given mass, known as the halo occupation number, (iV ga i)(m), and 
the second and third moments of the galaxy distribution ((N ga i(N ga \ — 1)} and (N ga i(N ga \ — l)(iV ga i — 2))). Following the method in 
Scoccimarro et al. (2001), the real-space halo model expressions for the 1-halo contributions to the galaxy 2PCF and 3PCF are given by 

£f£(r) = 2tt fdmn(m) ( N **W& - 1)) j dg j d6 > s 2 sin6 i Wm ( s ) Mm (| s + r |) ; 

J n gal Jo Jo 

<£{r,q,il>) = /dmn(m) <iVgal(iVgal T3 1)(iVgal ~ 2)) f "'da f d V f d6 s 2 sin 9 u m (s)u m (\s + r\)u m (\s + q\), (47) 

J n gal Jo JO J 

where n a denotes the mean number density of galaxies defined as 



dmn(m){N ga .i)(m). (48) 



In a similar manner, we can derive halo model predictions for the 2-halo contribution to the galaxy 2PCF, and the 2- and 3-halo contributions 
to the galaxy 3PCF. For simplicity we assume that the distribution of galaxies within a given halo follows the dark matter profile u m (r). 
Note that in the large-scale limit, the galaxy bias parameter is defined as 

bg = _L [dm n(m)b(m) <iYgal)(m) . (49) 
n a J m 

The halo occupation number has been used in the literature to explain galaxy clustering properties (Jing, Mo & Borner 1998; Seljak 
2000; Peacock & Smith 2000; Scoccimarro et al. 2001; Sheth et al. 2001; Berlind & Weinberg 2002; Moustakas & Somerville 2001; Scranton 
2002; Cooray 2002). In this paper, we employ the model of Scranton (2002), which approximately reproduces the results for (iV ga i) (m) from 
semi-analytic models in Kauffmann et al. (1999). The model gives {N ga .i) for red and blue galaxies 2 

f (-T-Y R e-^Ro/m) 1 ' 2 for red galaxies, 

(AT gal )(m) = M m *' A V S (50) 

\ (^) + Ae-M l °z™-™B 3 ) for blue galaxies 

where the parameters jr = 1.1, tur = 1.8 x 1O 13 /i _1 M and turo = 4.9 x 10 12 ft _1 M Q for red galaxies, while = 0.93, ms = 
2.34 x 10 13 /i _1 Mq, A = 0.65, Aq = 6.6 and ms 3 = 11.73 for blue galaxies. In the following, we employ a lower mass cutoff of 
m > 10 11 hT 1 Mq for the calculations, since in small halos effects such as supernova winds can blow away the gas from halos and thus 
suppresses further star formation. The model above yields bias parameters b giIcd = 1.51 and & 9i biuo = 0.82 for red and blue galaxies, 
respectively. 

Figure 16 plots the mean number of galaxies, (N gal }(m), as a function of parent halo mass. The figure shows that (iV ga i) < 1 at 
m 10 13 ft. -1 M© for both types of galaxies; a significant fraction of low mass halos are thus likely to contain at most one galaxy. For such 
low mass halos, we need to carefully model the higher order moments of the galaxy distribution within them, as discussed below. It is also 
apparent that (N gal ) for blue galaxies has a bump feature for halos of m £ 10 12 hT x M Q . The bump is due to the prescriptions used in the 
semi-analytical model. Although stars are formed from the cold gas in halos within a dynamical time scale, it is assumed that the gas cooling 
in halos with circular velocity V c > 350 km s _1 does not form visible stars so that the output galaxy catalog can fit the observed Tully-Fisher 
relation (see Kauffmann et al. 1999 for details). This sharp cutoff of the star formation prevents further formation of galaxies in halos with 
V c > 350 km s _1 unless the halos experience merging, which is taken as the trigger to induce the starburst. Note that the circular velocity 



2 Both samples are brighter than M v = — 17.7 + 51og/i in the simulations and the red and blue galaxies are labelled as being redder or bluer than 
B- I = 1.8. 
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Figure 16. Mean number of red and blue galaxies as a function of parent halo mass using the forms in equation (50). The thin dashed curve shows the model 
for blue galaxies without the bump feature, which we will use to demonstrate the effect on the halo model prediction for the 3PCF in Figure 18. 



V c = 350km s -1 roughly corresponds to halo of m ~ lO 13 A/0 at the present epoch; the velocity is larger for the same mass halo at earlier 
epochs. In fact, we will show that the bump feature drastically affects the configuration dependence and amplitude of the reduced 3PCF of 
blue galaxies. 

As mentioned above, we also need the second and third moments of the galaxy distribution within parent halos, {N ga x(N sa ,\ — 1)} 
and (N ga x(N ga i — l)(AT ga i — 2)}. We follow Scranton (2002), who expressed them as (Ar ga i(iV ga i — 1)) ~ o? g {m) (iV ga i) 2 and 
(iV ga i(A^ g ai - l)(iVgai - 2)) ^ a«(m)(iV gal } 3 , where 

, x J 1 m>lO 13 ^- 1 M , 

" sW _ \ In( > /m/lO 11 /i- 1 M ) m < W^h^Me. ( ' 

This model implies that galaxies follow a Poisson distribution for massive halos m Si 10 13 h -1 Mq and a sub-Poisson distribution for smaller 
halos. Other treatments of the sub-Poisson distribution have been proposed by Sheth et al. (2001) and Scoccimarro et al. (20001). 

We expect that the models of the 2nd and 3rd moments of N sa .\(m) will affect Q ga i at small scales which are dominated by the 1- 
halo term. Large scales Si 5 /i _1 Mpc, in the quasi-linear regime, are dominated by the 3-halo term. The behavior of Q ga i on the large 
scales can be roughly expressed in terms of the bias parameter (49) as Q ga i ~ Qmnss/bg. This follows from Q ga i ~ Cgai/(3C| a i) an d that: 
Cgai ~ bgCmass and £ ga i ~ 6g^ mass . As a result, the anti-biasing and bias of blue and red galaxies should lead to larger or smaller amplitudes 
of Q ga i compared to Q maS s, respectively. 

Figure 17 shows Q ga i for the galaxy 3PCF with triangle configurations parameterized by r, q and ip as in Figure 5. For red galaxies, the 
halo model combined with the halo occupation number leads to a weaker configuration dependence and smaller amplitude for Q ga i than the 
mass 3PCF, as suggested from the measurements in Jing & Borner (1998). On the other hand, the 3PCF of blue galaxies displays complex 
features. This is mainly due to the bump feature in (JV ga i)(m) at 10 11 < m S, 10 13 h~ 1 M ( v ) , as shown below. Hence, a detailed knowledge 
of the configuration dependence of Q ga i could be used to quantitatively constrain the occupation number, or more generally, the galaxy 
formation scenario as a function of observed galaxy properties (morphology, luminosity and so on) within a halo of given mass. 

The sensitivity of the input model of the halo occupation number to the galaxy 3PCF can be more clearly understood from Figure 18, 
which shows the Q ga i parameters for equilateral triangles. The 1-, 2- and 3-halo contributions to the galaxy 3PCF are explicitly shown. This 
figure may be compared with the result for the galaxy bispectrum shown in Figure 6 of Scoccimarro et al. (2001). It is again apparent that the 
amplitude of Q ga i for red galaxies is suppressed compared to the mass 3PCF, while the result for blue galaxies exhibits complex features. 
Figure 19 shows halo model predictions for the 2PCF and 3PCF of mass and blue and red galaxies separately. It is worth noting that the 2PCF 
and 3PCF of blue galaxies are well approximated by power laws at scales r Si 0.2 /i _1 Mpc, unlike the case for the mass or for red galaxies. 
This is due to the fact that galaxy formation is inefficient in massive halos and thus suppresses the 1-halo term for blue galaxies compared to 
the mass or to red galaxies which are formed by mergers. Also, Q ga i for blue galaxies has contributions from the 2- and 3-halo terms over a 
wider range of scales compared to that of red galaxies. This is mainly due to the bump feature in {N ga i) at low mass scales in Figure 16, and 
partly due to the shallower mass slope in (N sa x) oc m? . 

The lower panel of Figure 18 explores the origin of the complex features in Q ga i for blue galaxies at 0.1 Sj r 2 h~ 1 Mpc. The dot- 
dashed curve shows how the halo model prediction for Q ga i for blue galaxies changes if we suppress the bump feature in (iVg a i) at low 
mass scales by setting A = in equation (50). The model (AT ga i) used is shown by the thin dashed curve in Figure 16. One can see that 
the complex features in Q ga i disappear, except for the problematic bump feature at r ~ 2Mpc (see the discussion around Figure 8), and its 
amplitude becomes smaller than the mass 3PCF because of the change of the bias parameter to b g = 1.24 from 6 g = 0.82, which results in 
behavior similar to the red galaxies. Thus, the input model of {N ga .\) has a drastic effect on Q ga i. We also investigate how the halo model 
predictions are affected by altering the sub-Poisson distribution for the galaxy third-order moment, which is relevant for the range {N sa \) < 1 
or equivalently m Si 10 13 H~ 1 Mq. Scoccimarro et al. (2001) proposed the binomial distribution, where the third-order moments is given by 
(-Ngai(-Ngai — l)(AT ga ] — 2)) = ag(2a 2 — l)(iV ga i} 3 . Note that this model leads to negative values for the third moment for q 9 1. The 
thin dotted curve in Figure 18 shows the result of using this prescription — it changes Q ga i only for scales r Si 0.2 ft _1 Mpc. 
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Figure 17. Halo model predictions for the galaxy 3PCF as a function of triangle shapes parameterized by r, q and ip (see Figure 4). The solid, dashed and 
dotted curves show the results for the mass and for blue and red galaxies, respectively. We have employed the model of equation (50) for the the halo occupation 
number of red and blue galaxies. 



4.3 The 3-point correlation function of weak lensing fields 

Next we consider the halo model predictions for the angular 3PCF of the weak lensing convergence field. It should be noted that the following 
method can be easily extended to the predict the angular 3PCF of the galaxy distribution, as discussed in §3. 
The 1-halo term for the 3PCF of k is given by the following 4-dimensional integral: 

Z lh (r,q,4>) = / dxd 2 A (x) dmn(m;x) / d6 / dip 6 K m {6)n m {\0 + r \)K m (\& + q\), (52) 
Jo J Jo Jo 

where \0 + r | = (8 2 + r 2 + 26r cos %p) 1/2 and \ + q\ = (8 2 + q 2 + 26>gcos(v? - ip)) 1/2 - For the NFW profile, K m is analytically given 
by equations (26) and (28). Note that r and q denote angular scales, although we use the same notations as for the 3D case for simplicity. In 
analogy with equation (40), the Q parameter for the angular 3PCF is defined as 

Qk r,q,t/> Wr .(r)w K {q) + w K (r)w K Qr - q\) + w K (q)w K (\r - q\) ' 

where Z — Zih + Zih + Z^h- Note that the 2- and 3-halo terms, Z2h and Zzh, are computed from the approximations of equations (A6) 
and (A7), respectively. 

Figure 20 shows halo model predictions for Q K with triangle configurations parameterized as in Figure 5. The halos are taken to have 
NFW profiles. The hierarchical ansatz, Q K = const., does not hold rigorously over the scales we have considered, as in the 3D case. 
However, the configuration dependence of Q K for r = 0.1' is much weaker than for r — 1' and 10'. The strong configuration dependence in 
the two panels on the right can be understood as follows. First, the results in Figure 5 show that the Q parameter for the 3D mass distribution 
has a stronger configuration dependence for r £ 1 h~ Mpc than on strongly non-linear scales. The weak lensing convergence is a weighted 
projection of the mass distribution. As a result of the lensing projection, a range of 3-dimensional length scales contributes to Q K at a given 
0, which is therefore more sensitive to the 2- and 3-halo terms than the 3D case. This is explicitly shown in Figure 21, which plots the 1- 
and 2- and 3-halo contributions to Q K for the triangle shapes in Figure 20. For the smallest scale, r = 0.1', the 1-halo term indeed yields 
the dominant contribution to Q K , but the 2-halo term still has a non-negligible contribution (~ 10%). For r = 1', the 2-halo term becomes 
important ( £ 20% contribution). The 3-halo term becomes dominant at for r £ 10', at which Q K displays a characteristic, oscillatory feature 
as predicted by perturbation theory. 

In Figure 22 we show Q K for equilateral triangles as a function of side length. One can readily see that projection effects lead to greater 
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Figure 18. The upper panel shows the reduced 3PCF of red (dashed curves) and blue (dotted curves) galaxies as a function of the side length of equilateral 
triangles, as in Figure 7. The thin curves show the 1-, 2- and 3-halo contributions separately. For comparison, the solid curve shows the result for the mass 3PCF. 
Note that the bias parameters are b g rcd = 1.51 and b g blue = 0-82 for the red and blue galaxies, respectively. The lower panel explores the dependence 
of Qg a i for blue galaxies on the model ingredients. The dot-dashed curve shows the result for the Q ga i of blue galaxies, if we suppress the bump feature in 
(A r ga i)(m) at 10 11 < m 10 13 /i _1 A/q for the calculation. The thin dotted curve shows how the model of Scoccimarro et al. (2001) for the sub-Poisson 
distribution of blue galaxies changes the prediction at r ^ 0.2 ft~ 1 Mpc. 
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Figure 19. The 3PCF (left) and 2PCF (right) of mass, red and blue galaxies. 
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Figure 20. Halo model predictions for the reduced Q parameter for the weak lensing convergence, with triangle configurations parameterized as in Figure 5. 
The configuration dependence of Q changes as one shifts from the strongly nonlinear regime shown in the left panel to the quasi-linear regime shown in the 
right panel. 
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Figure 21. The 1-, 2- and 3-halo contributions to the Q parameter are shown for the convergence, with triangle shapes parameterized as in Figure 20. 



contributions from the 2- and 3-halo terms, in contrast to the 3D case shown in Figure 7. One caveat that should be noted, as discussed 
for Figure 7 and 8, is that the halo model might overestimate Q at the transition scale between the quasi-linear and non-linear regimes. 
This would be reflected in Q K over the angular scales 0.1' ^ r ^ 20', where the 2-halo term is relevant. This uncertainty needs further 
investigation by comparison with ray-tracing simulation results. The figure also shows that, on very small scales ^ 1', the 1-halo term of Q 
decreases slightly with decreasing angular scale. The small scale slope of the 1-halo term is determined by the combined effects of the NFW 
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r [arcmin] 

Figure 22. The Q parameter for the convergence for equilateral triangles as a function of the side length r (in arcminutes), as in Figure 7 and 9. 

profile, the mass function and the halo concentration in the halo model calculation. For lensing, the small scale slope of the 2PCF and 3PCF 
is not transparent due to projection effects. 

In the lower panel of Figure 22, we plot Q K for different cosmological models as shown in Figure 9. Q K is sensitive to cosmological 
parameters, in particular to Q m o parameter, in analogy with the sensitivity of the skewness parameters of weak lensing to Q mu (see e.g. 
Bernardeau et al. 2003; TJ02). Hence, it is expected that measuring Q K can be used to break the degeneracies in the determinations fl m o and 
as so far from measurements of the shear 2PCF (e.g., Van Waerbeke et al. 2001). 

4.4 Performance of our approximation for calculating the 1-halo term of the 3PCF 

In Appendix B we present a useful approximation for calculating the one halo contribution to the mass 3PCF, which uses the Fourier-space 
halo model combined with the approximation developed in TJ02 for calculating the skewness and kurtosis parameters. Figure 23 demonstrates 
the performance of our approximation by comparing the results with the exact values computed from the real-space halo model. We have used 
the approximation for the 3PCF calculation in the numerator of Q, and the 2PCF calculation in the denominator includes full contributions 
from the 1- and 2-halo terms. We consider some of the triangle shapes shown in Figure 5. The left and right panels show the results for 
r = 0.02 and 0.1 h~ 1 Mpc. In each panel, the solid, dashed and dotted curves are for q/r = 1, 3 and 5, respectively. Note that all the 
triangles considered are in the strongly non-linear regime since we wish to check the validity of our approximation for the 1-halo term. The 
thin curves are the exact predictions from the real-space halo model computed from equation (39), while the thick curves denote the results 
of the approximation (B4). The comparison demonstrates that the approximation works remarkably well with an accuracy of 5% for 
q/r = 1 and 3, and ^ 10% even for elongated triangles with q/r — 5. The three dot-dashed curves show the results of each term in the 
approximation (B4) for q/r — 3. Interestingly, although each term does not work well, the approximation obtained by summing them is 
close to the exact value. We have confirmed that the approximation of equation (B5) for the 1-halo term of the weak lensing 3PCF also works 
to about the same accuracy as the 3D case shown here. 



5 DISCUSSION AND CONCLUSION 

In this paper, we have used the halo clustering model to compute the 3-point correlation function (3PCF) of cosmological fields. We have 
shown results for the three-dimensional mass and galaxy distributions and the two-dimensional weak lensing convergence field. The halo 
model enables us to separately consider the contributions to the 3PCF arising from triplets in a single halo or two or three different halos. 
Thus we can understand how gravitational clustering transitions from the quasi-linear regime to the strongly non-linear regime as one goes 
to smaller scales. We found that the single halo contribution, which is dominant on small scales, can be computed using the real-space 
formulation of the halo model far more easily than the Fourier-space approach used in the literature. We also developed approximations 
for computing the 2- and 3-halo contributions to the 3PCF. Since measuring the real space 3PCF on small scales is likely to be easier than 
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Figure 23. The performance of our approximation (B4) for calculating the 1-halo contribution to the mass 3PCF. The results are shown for the triangle 
configurations as in Figure 5. We consider the non-linear scales r = 0.02 (left panel) and 0.1 h _1 Mpc (right). The thin curves are the exact halo model 
predictions, while the thick curves show the results of the approximation of equation (B4). The three dot-dashed curves show each term in the approximation 
(B4) for q/r = 3. This figure verifies that our approximation accurately describes the 3PCF, even though each term in equation (B4) is far from the exact 
value. 



the bispectrum, our model predictions will allow for the extraction of cosmological information from forthcoming galaxy surveys and weak 
lensing surveys. 

We obtained the following results by applying our halo model to the concordance CDM model (ACDM). The 3PCF on small scales 
r H 1 ft~ 1 Mpc is dominated by the 1-halo term. Hence it probes properties of massive halos and can be used to constrain halo profiles as 
discussed below. The quasi-linear 3PCF predicted by perturbation theory can be reproduced by the 3-halo contribution for r ^ 3 fo _1 Mpc. 
For plausible halo model parameters, the hierarchical ansatz for the reduced 3PCF parameter, Q = constant, does not hold over the range 
of scales we have considered. However the Q parameter does have a weaker dependence on triangle configurations in the non-linear regime 
than in the quasi-linear regime (see Figures 5 and 10). These results qualitatively verify the results in Ma & Fry (2000b,c). Ma & Fry also 
pointed out that the halo model for plausible model parameters violates the stable clustering hypothesis which has been widely used to 
develop analytical prediction of the non-linear gravitational clustering, as done in the popular PD fitting formula. In fact, very recently Smith 
et al. (2002) showed a weak violation of the stable clustering hypothesis using high-resolution TV-body simulations. However, our halo model 
predictions for the 3PCF do not match the asymptotic shape proposed in Ma & Fry (2000b) because of the assumptions employed in the 
asymptotic formula. Therefore, it it merits further study to carefully investigate how small-scale gravitational clustering can be described by 
the halo model ingredients. Such a study can be carried out with the halo model methods developed in this paper. 

We also found that the non-linear 3PCF is most sensitive to the inner slope of the halo profile and to the halo concentration parameter 
for given cosmological parameters (see Figures 12, 13 and 14). Combinations of the inner slope and the concentration can be adjusted so 
that the halo model matches the 2PCF result (see Figure 2). However only one combination can match both the 2- and 3PCFs. Hence in 
combination with the 2PCF, the non-linear 3PCF could be used to constrain the properties of dark matter halos. We suggest that the use of 
higher order correlations is a useful way of measuring the properties of massive halos. For example, while lensing surveys of clusters can be 
used to measure cluster halo profiles directly, this involves identifying cluster centers and assigning masses to clusters to measure averaged 
profiles. In contrast, by measuring the 2- and 3PCFs, parameters of the halo mass function and profile can be fitted for. While this approach 
is more challenging computationally, and requires some theoretical assumptions, it treats the data more objectively. 

Our halo model predictions match the simulation results of Jing & Borner (1998) in the non-linear regime as well as the quasi-linear 
regime, as shown in Figure 10. However, the halo model seems to be less accurate on the intermediate transition scale ~ 1 /i -1 Mpc. Figure 7 
shows that the predicted Q parameter for equilateral triangles has a bump at this scale, which corresponds to the bump found for the reduced 
bispectrum at k ~ 1 h Mpc -1 by Scoccimarro et al. (2001). This is unlikely to reflect real properties of dark matter clustering. Rather, 
we argued that the bump feature is sensitive to the sharp cutoff of halos at the virial radius and the exclusion effect for the 2- and 3-halo 
terms. We explored alternative prescriptions as shown in Figure 8. To resolve this, detailed calculations and comparisons of the halo model 
predictions with simulation results will be needed. 

We extended the halo model to predict the 3PCF of the galaxy distribution. Once we model how galaxies populate a halo of given 
mass, the halo occupation number, we can straightforwardly predict the galaxy 3PCF based on the halo model developed here. For the halo 
occupation number for red galaxies we have used, the galaxy 3PCF has a smaller amplitude and weaker dependence on triangle configuration 
compared to the mass 3PCF (Figures 17 and 18). This trend is indeed consistent with the actual measurements in Jing & Borner (1998). 
On the other hand, the 3PCF of blue galaxies displays complex features reflecting the input model of the halo occupation number. Thus, the 
galaxy 3PCF can be used to constrain the galaxy formation scenario as a function of host halo properties and galaxy type. Further work is 
needed to model the expected properties of galaxies by type for specific surveys and compute the resulting 3PCF. 

As an example of the angular 3PCF, we have computed the 3PCF of the weak lensing convergence field. In particular, we employed 
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the real-space halo model to compute the single halo contribution, as for the 3D case, which enabled us to compute the 1-halo term by a 
4-dimensional integration. We verified that the real-space halo model is equivalent to the Fourier-space model well studied in the literature 
(see Figure 3). We also developed approximations for calculating the 2- and 3-halo terms. Because of projection effects, the 2- and 3-halo 
terms contribute to the 3PCF over a wider range of scales for the angular 3PCF compared to the 3D 3PCF. The resulting 3PCF does not obey 
the hierarchical ansatz over the angular scales we have considered. The lensing 3PCF is sensitive to cosmological parameters, in particular 
O m() . Comparing measurements with model predictions of the 3PCF can be a useful tool to break degeneracies between the power spectrum 
and Q m o. We intend to compute 3-point functions of the shear field following Schneider & Lombardi (2002) and Zaldarriaga & Scoccimarro 
(2002) as these are easier to measure from the data. In Takada & Jain (2003), we presented a brief investigation of the 3PCFs of shear fields, 
where we found good agreement between the halo model prediction and the results measured from the ray-tracing simulations. It is also 
expected that the n-point correlations of weak lensing on sub-arcminute scales contain a wealth of information on properties of massive 
halos, beyond their dependence on cosmological parameters. The real-space halo model we have developed will be a useful analytical tool 
for such calculations. 

The halo model presented in this paper allows for several interesting applications. First, the model can be extended to investigate the 
effect of triaxial halo shapes on the 3PCF, since it is the lowest-order statistical probe of non-sphericity. So far halo model applications assume 
a single, spherically symmetric profile. Recently, based on high-resolution simulations, Jing & Suto (2002) showed that halos are more 
accurately described by triaxial halo profiles than the spherically symmetric NFW profile. They claimed that the axis ratios typically have 
vales of ~ 0.6 for the smallest and largest axis. The non-sphericity of the halo profile could lead to a characteristic configuration dependence 
of the 3PCF. This effect should also affect the interpretation of cosmic shear measurements (Jing 2002). Likewise, the formulation developed 
in Sheth & Jain (2002) can be used to discuss the effects of substructure within halos on the 2PCF and 3PCF. 

The real-space halo model can be directly used to compute the covariance of the 2PCF. As discussed in Cooray & Hu (2001b), the mass 
distribution on small scales displays pronounced non-Gaussian features induced by non-linear gravitational clustering. Hence one needs to 
take into account the non-Gaussian errors arising from the 4-point correlation function (4PCF). On the scales of interest the single halo 
contribution should dominate the error. The real-space halo model allows us to analytically compute the error contribution arising from the 
4PCF with no additional computational effort than for the 3PCF (see §2.1). 

In Appendix B, we constructed an approximation for calculating the 1-halo term of the 3PCF based on the Fourier-space halo model. 
We showed that the the approximation accurately describes the amplitude and configuration dependence of the 3PCF (see Figure 23). We will 
employ this approximation to perform an analytical study of the pairwise peculiar velocity dispersion (PVD) within the BBGKY hierarchy 
formalism (Peebles 1980). The PVD can be measured through the redshift distortions inferred from galaxy surveys (e.g., Zehavi et al. 2002). 
The BBGKY picture tells us that the PVD arises mainly from the 3PCF on small scales. However, there has been no comprehensive analytical 
model to describe the non-linear PVD. This is because of the complex form of the BBGKY hierarchy equations. Peebles (1980) (see Mo, 
Jing & Borner 1997; Jing et al. 1998 for a detailed study) assumed the hierarchical form for the non-linear 3PCF, although it turns out to be 
violated for CDM models. The isothermal assumption for the velocity distribution within a given halo was employed to analytically obtain 
the PVD based on the halo model (Sheth et al. 2001). It is important to clarify whether the BBGKY hierarchy leads to a self-consistent PVD 
for the CDM model. Moreover, we can easily combine the halo model prediction with models of the halo occupation number of galaxies to 
predict the PVD for different galaxy types, following the approach in §4.2.3. 

We would like to thank R. Sheth, R. Scoccimarro, I. Szapudi, A. Taruya and L. Hui for several valuable discussions. We thank Y. P. 
Jing for kindly providing us with his simulation data results. Helpful comments from the referee, Chung-Pei Ma, led to improvements in the 
paper. This work is supported by NASA grants NAG5-10923, NAG5-10924 and a Keck foundation grant. 



APPENDIX A: APPROXIMATION FOR THE 2- AND 3-HALO TERMS OF THE 3PCF 

In this appendix, we present approximations used for the predictions of the 2- and 3-halo contributions to the 3PCF of the density field. 
The approximations are based on the Fourier-space halo model, combined with an approximation for the configuration dependence of the 
bispectrum from Scoccimarro et al. (2001; see also TJ02). 



Al 3D 3-point correlation function 



Following the Fourier-space halo model, as discussed in §2.2, the 2-halo term for the mass 3PCF can be expressed from equations (11) and 
(13) as 



(2h(ri,r 2 ,r 3 ) 
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(2tt) 3 
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d 3 k 2 



(277)3 



' (k 2 )P (fe) exp[ife 2 • r 32 ] 



+ eye. 



(Al) 



where the bias parameter b(m) is given by equation (20). As discussed in Scoccimarro et al. (2001), the physical meaning of the 2-halo term 
tells us that the main contribution to the first term on the r.h.s arises from modes with ki 3> k 2 , leading to the approximation ki 2 ~ fci. 
Replacing u m {ki 2 ) with u m (ki) allows us to perform analytically the angular integration in j d :i ki in the first term above. Using similar 
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procedures for the second and third terms on the r.h.s. yields the following approximation: 
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where we have used the formula fdflfe cxp[ifc ■ r] = jo (kr). The large square bracket is meant to show explicitly that the integrations inside 
each bracket can be done separately from the other. Thus one needs to perform only 2-dimensional integrations to get the 2-halo term. It was 
shown that a similar approximation for the skewness calculation turns out to be accurate (Scoccimarro et al. 2001; TJ02). Hence, we expect 
that this also holds for the 3PCF calculation. 

The 3-halo contribution to the 3PCF can be expressed as 
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We again employ the approximation £t m // (fci 2 ) ~ u m u (ki) to perform analytically the integration for the first term on the r.h.s., although 
there is no guarantee that the approximation is as accurate as for the 2-halo term. Following the procedure developed to derive the perturbation 
theory 3PCF (Fry 1984b; also see Jing & Borner 1997; Gaztanaga & Bernardeau 1998; Barriga & Gaztanaga 2002) we can derive the 
approximation for the 3-halo term: 
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We find that this approximation reproduces the perturbation theory result at scales £ 5 h~ 1 Mpc for the ACDM model (see Figure 7). Note 
that, for the actual predictions of the 3PCF shown in this paper, we have used the permutation (1) <-» (2) for terms such as E^E' 1 - 1 in the 
above equation, so that it satisfies statistical symmetry for permutations between n, r 2 and r 3 . 

By replacing (m/po) in equations (A2) and (A4) with the halo occupation number of galaxies {N gal )(m)/n S!l i, as discussed in §4.2.3, 
we obtain the corresponding approximations for computing the 2- and 3-halo contributions to the 3PCF of galaxies. 



A2 2D 3-point correlation function 



It is straightforward to extend the approximations for the 3-dimensional 3PCF to the 3PCF of weak lensing fields. As discussed for equa- 
tion (36), combining equation (A2) with Limber's equation yields the following approximation for the 2-halo term of the weak lensing 
convergence: 
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where ki — /;/eU(x)- Note that the redshift dependence of P L (k; x) is given by the linear growth factor D(z). 

Likewise, the approximation for the 3-halo term can be obtained by the following replacements in equation (A4) and the additional 
integration over the redshift: 
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where ^1323 = {ri3 ■ r 23 )/(ri3r23) and 
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APPENDIX B: APPROXIMATIONS FOR THE 1-HALO TERM OF THE 3PCF 

In this section, we present approximations for calculating the 1-halo contribution to the 3D or 2D 3PCF with the Fourier-space halo model. 
These are combined with the approximation in TJ02 for calculating the skewness parameter of the density field. 

Bl 3D 3-point correlation function 

Based on the Fourier-space halo approach, the 1-halo contribution to the 3PCF of the density field can be obtained from equations (11) and 
(6) as 

Cih{ri,r2,r 3 ) = jdmn(m)f^-j /jj ^~^s"™( fei ) ex P[* fc i • ri](27r) 3 5r>(fei23). (Bl) 

%— i 

The Dirac-delta function allows us to eliminate one of the three wavevectors fc;. For example, let us eliminate fc3, leading to the replacement 
k :i — — ki2 in the above equation. Then, the integrand has dependences on d 3 ki via u m (fci2) and exp[ifc; • r;]. TJ02 proposed that the 
angular dependence of u m (fci 2 ) can be approximated with the replacement: 

u m (k 12 ) — > u m (k 12 ), (B2) 

where ki2 = (kf + k\ — k\k 2 ) x ^ . TJ02 showed that this replacement works well for predictions of the skewness parameters of the 3D 
the 

integrations in j d 3 ki, to get 



density field and the weak lensing field 3 . If we assume that this holds for the 3PCF calculation, we can analytically perform the angular 



Cih(fi,r 2 ,r 3 ) « F c (ri 3 ,r 23 ) = Jdmn(m)\^-^j Jt^t^ k i klu m {ki)uj n (k2)u m {ki2)3o{kiriz)jo{k2r2z)- (B3) 

The above approximation reduces the 7-dimensional integral in equation (Bl) to a 3-dimensional one. However, the resulting 3PCF depends 
only on the two side lengths ri3 and r 2 3, and thus has no dependence on the interior angles. If we had started by eliminating k 2 or fei in 
equation (Bl), the resulting 3PCF would have dependences only on the two side lengths ri2 and r23 or ri2 and ri3, respectively. Thus, the 
approximate result for the 3PCF has different forms even for the same triangle shape. To resolve this, we propose the following symmetrized 
form as the approximation for calculating (i^: 

CiT( r i^2,r 3 ) = i [F c (ri2,r 13 ) +F f (ri2,r 2 3) +F f (r 13 ,r2 3 )] , (B4) 

The dependences of this form on the three parameters of the triangle shape, ri2, ^23 and r3i, are like the hierarchical form £ = 
<3[£(ri 2 )£(?*3i) + £( r i2)£(r23) + £(ri 3 )£(r22)]- The factor 1/3 is simply chosen so that it agrees with the third-order moment given in 
TJ02 for the limit r, q — > 0. As shown in §4.4, the approximation (B4) is accurate compared with the exact value computed from equation 
(6), even though each term on the r.h.s is far from the exact value (see Figure 23). 

B2 2D 3-point function 

In the spirit of equation (B4), the corresponding approximation for the projected 3PCF given by equation (35) can be constructed as: 

^(01,02,03) = i [F Z (e i2 , e i3 )+F Z (021, 023) +F Z (031, 9 3 2)] , (B5) 

with 

F z (0,<p) = £ S d X W :i (x)d^( X ) jdmn(m)(j^j J u m (ki)u m (k2)u m (ki2)J O (0li)Jo(yl2), (B6) 

where fc» = U/dA(x)- The 4-dimensional integration required to carry out the above approximation is the same level of complexity as the 
real-space halo model expression in equation (35). 

3 Note that, although one might imagine another possible replacement u m (fci2) — » Um{ki), it overestimates the true value by about ~ 50% for the 1-halo 
term, as shown in TJ02. 
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